library(MCMCpack)
library(coda)
library(ggplot2)

data<-read.csv("spanish_tc_votes_for_export.csv",header=T)

## Set up theta constraints for MCMCirt1d
## These must be done by column number
## Our variable names are a bit too long otherwise
## We set Garvia-Calvo to + because he was a figure under the Francoist regime
## We then set Perez-Tremps to - because he is often reported as left-leaning

my.theta.constraints<-list(V19="+",V14="-") 

iterations<-1.75e6
burnin<-.25e6

votes<-t(data[,8:ncol(data)])
rownames(votes)<-NULL

## Run the model
## This may take a while
set.seed(1982)
my.model<-MCMCirt1d(votes,
    theta.constraints=my.theta.constraints,
    mcmc=iterations,
    burnin=burnin,
    thin=250,verbose=iterations/3,
    drop.constant.items=FALSE,
    store.item=TRUE)

## Standardize our ideal points to have mean zero and SD 1

ideal.points<-grep("theta",colnames(my.model))
difficulties<-grep("beta",colnames(my.model))
locations<-grep("alpha",colnames(my.model))

normalize.model<-function(mat) {
	my.mat<-mat
	ref.mean<-mean(my.mat[ideal.points],na.rm=T)
	ref.sd<-sd(my.mat[ideal.points],na.rm=T)

	my.mat[ideal.points]<- (my.mat[ideal.points] - ref.mean)/ref.sd

	my.mat[locations]<- my.mat[locations] - ref.mean * my.mat[difficulties]

	my.mat[difficulties]<- my.mat[difficulties] * ref.sd

	my.mat
}

my.model<-as.mcmc(t(apply(my.model,1,normalize.model)))

## Calculate measures of fit, etc.,
source("diagnostics.r",echo=TRUE)

## Plot judges
results.df<-data.frame(holder$quantiles[grep("theta",rownames(holder$quantiles)),])
names(results.df)<-c("Lo","X5","Median","X95","Hi")
results.df$Judge<-colnames(data[,8:ncol(data)])

party<-read.csv("spanish_judges.csv",header=T)
results.df$Party<-party$Party[match(party$Judge,results.df$Judge)]
results.df$DisplayName<-party$DisplayName[match(party$Judge,results.df$Judge)]

results.df$Party<-factor(results.df$Party,ordered=TRUE)
results.df$Party<-reorder(results.df$Party,results.df$Median)
my.order<-order(results.df$Party,results.df$Median)
results.df$DisplayName<-factor(results.df$DisplayName,levels=results.df$DisplayName[my.order],ordered=TRUE)

ideal.plot<-ggplot(data=results.df,aes(x=DisplayName,y=Median,ymax=Hi,ymin=Lo,color=Party,shape=Party)) + geom_pointrange() + theme_bw() + scale_y_continuous("Position") + scale_x_discrete("Judge") + coord_flip()

pdf(file="spanish_ideal_points.pdf")
ideal.plot
dev.off()

write.csv(results.df,"spanish_ideal_points.csv",row.names=FALSE)

## Plot discriminating cases by case type
cases.df<-data.frame(
	as.vector(by(my.discriminating,data$Basis,mean)),
	levels(data$Basis),
	as.vector(table(data$Basis)))
names(cases.df)<-c("Proportion","Basis","Number")
cases.df$Basis<-as.character(cases.df$Basis)
cases.df$Basis<-paste(cases.df$Basis, " (",cases.df$Number,")",sep="")

pdf(file="spanish_case_types.pdf")
ggplot(cases.df,aes(x=Proportion,y=Basis)) + geom_point(size=2.5) + theme_bw(base_size=18) + scale_x_continuous("% of discriminating cases") + opts(axis.title.y=theme_blank())
dev.off()
save.image("spanish.RData")
